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Abstract 

Background. Inferring the structure of gene regulatory networks (GRN) from a col- 
lection of gene expression data has many potential applications, from the elucidation of 
complex biological processes to the identification of potential drug targets. It is however a 
notoriously difficult problem, for which the many existing methods reach limited accuracy. 

Results. In this paper, we formulate GRN inference as a sparse regression problem and 
investigate the performance of a popular feature selection method, least angle regression 
(LARS) combined with stability selection, for that purpose. We introduce a novel, robust 
and accurate scoring technique for stability selection, which improves the performance of 
feature selection with LARS. The resulting method, which we call TIGRESS (for Trustful 
Inference of Gene REgulation with Stability Selection), was ranked among the top GRN 
inference methods in the DREAM5 gene network reconstruction challenge. We investigate 
in depth the influence of the various parameters of the method, and show that a fine param- 
eter tuning can lead to significant improvements and state-of-the-art performance for GRN 
inference. 

Conclusions. TIGRESS reaches state-of-the-art performance on benchmark data. This 
study confirms the potential of feature selection techniques for GRN inference. 

Availability. Code and data are available on http:/ /cbio.ensmp.fr/^ahaury. Moreover, 
running TIGRESS online is possible on GenePattern: http://www.broadinstitute.org/ 
cancer /sof tware/genepattern/. 

1 Background 

In order to meet their needs and adapt to changing environments, cells have developed vari- 
ous mechanisms to regulate the production of the thousands of proteins they can synthesize. 
Among them, the regulation of gene expression by transcription factors (TF) is preponderant: 
by binding to the promoter regions of their target genes (TG), TF can activate or inhibit their 
expression. Deciphering and understanding TF-TG regulations has many potential far-reaching 
applications in biology and medicine, ranging from the in silico modelling and simulation of the 
gene regulatory network (GRN) to the identification of new potential drug targets. However, 
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while many TF-TG regulations have been experimentally characterized in model organisms, the 
systematic experimental characterization of the full GRN remains a daunting task due to the 
large number of potential regulations. 

The development of high-throughput methods, in particular DNA microarrays, to monitor 
gene expression on a genome-wide scale has promoted new strategies to elucidate GRN. By 
systematically assessing how gene expression varies in different experimental conditions, one can 
try to reverse engineer the TF-TG interactions responsible for the observed variations. Not 
surprisingly, many different approaches have been proposed in the last decade to solve this GRN 
reverse engineering problem from collections of gene expression data. When expression data 
are collected over time, for example, several methods have been proposed to construct dynamic 
models where TF-TG interactions dictate how the expression level of a TG at a given time 
allows to predict the expression levels of its TG in subsequent times [2, 21, 9, 1, 36, 33, 16, 8, 
5, 4]. When expression data are not limited to time series, many methods attempt to capture 
statistical association between the expression levels of TG and candidate TF using correlation 
or information-theoretic measures or mutual information [7, 26, 11] or take explicit advantage 
of perturbations in the experiments such as gene knock-downs [31]. The difficulty to separate 
direct from indirect regulations has been addressed with the formalism of Bayesian networks 
[14, 17, 13], or by formulating the GRN inference problem as a feature selection problem [19]. 
We refer to [27, 24] for detailed reviews and comparisons of existing methods. 

Recent benchmarks and challenges have highlighted the good performance of methods which 
formalize the GRN inference problem as a regression and feature selection problem, namely, 
identifying a small set of TF whose expression levels are sufficient to predict the expression level 
of each TG of interest [28]. This idea underlies the Bayesian network formalism [14], but is 
more directly addressed by GENIE3 [19], a method which uses random forests to identify TF 
whose expression levels are predictive for the expression level of each TG, and which is now 
recognized as state-of-the-art on several benchmarks [19, 24]. Feature selection with random 
forests remains however poorly understood theoretically, and one may wonder how other well- 
established statistical and machine learning techniques for feature selection would behave to 
solve the GRN inference problem. 

In this paper, we investigate the performance of a popular feature selection method, least 
angle regression (LARS) [ I ] combined with stability selection [3, 29], for GRN inference. LARS 
is a computationally efficient procedure for multivariate feature selection, closely related to 
Lasso regression [3 1]. We introduce a novel, robust and accurate scoring technique for stability 
selection, which improves the performance of feature selection with LARS. The resulting method, 
which we call TIGRESS (for Trustful Inference of Gene REgulation with Stability Selection), was 
ranked among the top GRN inference methods in the DREAM5 gene reconstruction challenge 
[23]. We furthermore investigate in depth the influence of the various parameters of the method, 
and show that a fine parameter tuning can lead to significant improvements and state-of-the-art 
performance for GRN inference. Overall this study confirms the potential of state-of-the-art 
feature selection techniques for GRN inference. 

2 Methods 

2.1 Problem formulation 

We consider a set of p genes Q = [l,p], including a subset 7~ C [l,p] of transcription factors, 
among which we wish to discover direct regulations of the form (t, g) for t £ T and g € Q. 
We do not try to infer self-regulation, meaning that for each target gene g £ Q we define the 
set of possible regulators as T g = T\{g} if g £ T is itself a transcription factor, and Tg = T 
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otherwise. The set of all candidate regulations is therefore £ = {(t,g),g E Q,t E T g }, and the 
GRN inference problem is to identify a subset of true regulations among £. 

For that purpose, we assume we have gene expression measurements for all genes Q in n 
experimental conditions. Although the nature of the experiments may vary and typically include 
knock-down or knock-out experiments and even replicates, for simplicity we do not exploit this 
information and only consider the nxp data matrix of expression levels X as input for the GRN 
inference problem. Each row of X corresponds to an experiment and each column to a gene. We 
assume that the expression data have been pre-processed for quality control and missing values 
imputation. 

In order to infer the regulatory network from the expression data X, we compute a score 
s : £ — > M. to assess the evidence that each candidate regulation is true, and then predict as true 
regulation the pairs (t,g) E £ for which the evidence s(t,g) is larger than a threshold S. We let 
5 as a user-controlled parameter, where larger 5 values correspond to less predicted regulations, 
and only focus on designing a significance score s(t,g) that leads to "good" prediction for some 
values of 8. In other words, we only focus on finding a good ranking of the candidate regulations 
£, by decreasing score, such that true regulations tend to be at the top of the list; we let the 
user control the level of false positive and false negative predictions he can accept. 

2.2 GRN inference with feature selection methods 

Many popular methods for GRN inference are based on such a score. For example, the correlation 
or mutual information between the expression levels of t and g along the different experiments is 
a popular way to score candidate regulations [7, 26, 11]. A drawback of such direct approaches is 
that it is then difficult to separate direct from indirect regulations. For example, if t\ regulates £2 
which itself regulates g, then the correlation or mutual information between t\ and g is likely to 
be large, although (t\,g) is not a direct regulation. Similarly, if t\ regulates both £2 arid g, then 
ti and g will probably be very correlated, even if there is no direct regulation between them. In 
order to overcome this problem, a possible strategy is to post-process the predicted regulations 
and try to remove regulations likely to be indirect because they are already explained by other 
regulations [26]. Another strategy is, given a target gene g E Q, to jointly estimate the scores 
s(t,g) for all candidate regulators t E T g simultaneously, with a method able to capture the 
fact that a large score for a candidate regulation (t, g) is not needed if the apparent correlation 
between t and g is already explained by other, more likely regulations. 

Mathematically, the latter strategy is closely related to the problem of feature selection 
in statistics, as already observed and exploited by several authors [28, 19]. More specifically, 
for each target gene g E G, we consider the regression problem where we wish to predict the 
expression level of g from the expression level of its candidate regulators t € Tg. 

Xg = fg(X Tg )+e, (1) 

where X{ represents the expression level of the z-th gene across different experiments (modelled 
as a random variable), Xj- = {Xt , t E Tg} is the set of expression levels of the candidate 
transcription factors for gene g, and e is some noise. Any linear or nonlinear statistical method 
for regression can potentially be used to infer f g from the observed expression data. However, 
we are not directly interested in the regression function f g , but instead in the identification of 
a small set of transcription factors which are sufficient to provide a good model for X g . We 
therefore need a score s g (t) for each candidate transcription factor t E T g to assess how likely it 
is to be involved in the regression model f g . For example, if we model f g as a linear function 

f 9 (Xr a ) = Y,^t, (2) 
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then the score s g (t) should typically assess the probability that j3t g is non-zero [28]. More general 
models are possible, for example [19] model f g with a random forest [6] and score a predictor 
s g (t) with a variable importance measure specific to this model. Once a score s g (t) is chosen to 
assess the significance of each transcription factor in the target-gene-specific regression model 
(1), we can combine them across all target genes by defining the score of a candidate regulation 
(t,g) £ £ as s(t,g) = s g (t), and rank all candidate regulations by decreasing score for GRN 
inference. 

2.3 Feature selection with LARS and stability selection 

We now propose a new scoring function s g (t) to assess the significance of a transcription factor 
t £ T g in the regression model (1). Our starting point to define the scoring function is the 
LARS method for feature selection in regression [10]. LARS models the regression function (1) 
linearly, i.e. it models the expression of a target gene as a linear combination of the expression 
of its transcription factors, as in (2). Starting from a constant model where no TF is used, it 
iteratively adds TF in the model to refine the prediction of X g . Contrary to classical forward 
stepwise feature selection [35, 18], LARS does not fully re-optimize the fitted model when a 
new TF is added to the model, but only refines it partially. This results in a statistically sound 
procedure for feature selection, akin to forward stage- wise linear regression and the Lasso [34, 18], 
and a very efficient computational procedure. In practice, after L steps of the LARS iteration, 
we obtain a ranked list of L TF selected for their ability to predict the expression of the target 
gene of interest. Efficient implementations of LARS exist in various programming languages 
including R (lars package, [10]) and MATLAB (SPAM toolbox, [22]). Since the selection of TF 
is iterative, LARS has the potential to disregard indirect regulations. 

The direct use of LARS to score candidate regulations has, however, two shortcomings. 
First, LARS can be very sensitive and unstable in terms of selected features when there exist 
high correlations between different explanatory variables. Second, it only provides a ranking of 
the TF, for each TG of interest, but does not provide a score s g (t) to quantify the evidence that 
a TF t regulates a target gene g. Since we want to aggregate the predicted regulations across 
all target genes to obtain a global ranking of all candidate regulations, we need such a score. 

To overcome both issues, we do not directly score candidate regulations with the LARS, but 
instead perform a procedure known as stability selection [29] on top of LARS. The general idea 
of stability selection is to run a feature selection method many times on randomly perturbed 
data, and score each feature by the number of times it was selected. It was shown that stability 
selection can reduce the sensitivity of LARS and Lasso to correlated features, and improve 
their ability to select correct features [3, 29]. In addition, it provides a score for each feature, 
which can then be aggregated over different regression problems, i.e. different target genes in 
our case. More precisely, to score the candidate target genes t £ T g of a given target gene g 
using LARS with stability selection, we fix a (large) number of iterations R, and repeat R/2 
times the following iterations: we randomly split the experiments into two halves of equal or 
approximately equal size, we multiply the expression levels of the candidate transcription factors 
in T g on each microarray by a random number uniformly sampled on the interval [a, 1] for some 
< a < 1, and we run the LARS method for L > steps on the two resulting reduced and 
reweighed expression matrices. We therefore perform a total of R LARS runs on randomly 
modified expression matrices. For each run, the result of LARS after L steps is a ranked list of 
L TF. After the R runs, we record for each g £ Q, t £ T g and / £ [1, L] the frequency F(g,t, I) 
with which the TF t was selected by the LARS in the top I features to predict the expression 
of gene g. We divide the frequency by R to obtain a final score between and 1, 1 meaning 
that t is always selected by LARS in the top I features to predict the expression level of g, and 
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that is is never selected. Figure 1 displays graphically these frequencies, for a given gene g 
fixed, all candidate TF in Tg, and I = 1, . . . , 15. When I increases, the frequency F(g, t, I) for 
fixed g and t is non-decreasing because the LARS method selects increasing sets of TF at each 
step. In addition, since the total number of TF selected after I LARS steps is always equal to I, 
taking the average over the R LARS runs leads to the equality XlteT ^G?) ^ = ^ f° r an y g ene 
g and LARS step I. 




Number L of LARS steps 

Figure 1: Illustration of the stability selection frequency F(g, t, I) for a fixed target gene g. Each 
curve represents a TF t S Tg, and the horizontal axis represents the number I of LARS steps. 
F(g, t, I) is the frequency with which t is selected in the first I LARS steps to predict g, when the 
expression matrix is randomly perturbed by selecting only a limited number of experiments and 
randomly weighting each expression array. For example, the TF corresponding to the highest 
curve was selected 57% of the time at the first LARS step, and 81% of the time in the first two 
LARS steps. 



Once the frequency table F(g, t, I) is computed for I = 1, . . . , L, we need to convert it into a 
unique score s(t,g) for each candidate pair (t,g). The original stability selection score [3, 29] is 
simply defined as the frequency of selection in the top L variables, i.e., 

Soriginal (*, s) = F(g,t,L) . (3) 

As suggested by Figure 1, this score may be very sensitive to the choice of L. In particular, if 
L is too small, many TF may have zero score (because there are never selected in the top L 
TFs), but when L is too large, several TF may have the same score 1 because they are always 
selected in the top L TFs. To alleviate this possible difficulty, we propose as an alternative score 
to measure the area under each curve up to L steps, i.e. to consider the following area score: 



(t,g) = jJ2F(g,t,l). (4) 



L 

i=i 

The difference between s or i g i na i(t, g) and s area {t,g) becomes clear if we consider the rank of t in 
the list produced by LARS in one run as a random variable Ht (with Ht = 1 meaning that t is 
ranked first by LARS). F(g,t,l) is then, by definition, the empirical probability P(Ht < I) that 
Ht is not larger than I. The original score has therefore an obvious interpretation as P{Ht < Li), 
which we can rewrite as: 



Soriginal\t,g) = F [<1> original (Ht)] with <poriginal(h) 



1 if h < L , 
otherwise. 



5 



Interestingly a small computation shows that the area score has a similar probabilistic interpre- 
tation: 

L 

1=1 

*£P(H t <l) 

1=1 

1=1 h=l 
L 

Y^{L + 1 - h)P(H t = h) 

h=l 

E[(j) area {Ht)\ , 

(L + l-h if h<L, 
1 otherwise. 

In other words, both the original and the area scores can be expressed as E[(j)(Hf)], although 
with a different function 0. While the original score only assesses how often a feature ranks in 
the top L, the area score additionally takes into account the value of the rank, with features 
more rewarded if they are not only in the top L but also frequently with a small rank among 
the top L. Since s area integrates the frequency information over the full LARS path up to L 
steps, it should be less sensitive than s or i g i na i to the precise choice of L, and should allow to 
investigate larger values of L without saturation effects when several curves hit the maximal 
frequency of 1. We note that other scores of the form E [cfi(Ht)} for non-increasing function cj) 
could be investigated as well. 

2.4 Parameters of TIGRESS 

In summary, the full procedure for scoring all candidate edges in £, which we call TIGRESS, 
splits the GRN inference problem into p independent regression problems taking each target 
gene g £ Q in turn, and scores each candidate regulation (t,g) for a candidate TF t £ T g with 
the original (3) or area (4) stability score applied to LARS feature selection. In addition to the 
choice of the scoring method (original or area), the parameters of TIGRESS are (i) the number 
of runs R performed in stability selection to compute the scores, (ii) the number of LARS steps 
L, and (iii) the parameter a £ [0, 1] which controls the random re-weighting of each expression 
array in each stability selection run. Apart from R that should be taken as large as possible 
to ensure that frequencies are correctly estimated, and is only limited by the computational 
time we can afford to run TIGRESS, the influence of a and L on the final performance of the 
method are not obvious. Taking a = 1 means that no weight randomization is performed on the 
different expression arrays, while a = leads to maximal randomization. [29] advocate that a 
value between 0.2 and 0.8 is often a good choice. Regarding the choice of L, [ )] mentions that 
it has usually little influence on the result, but as discussed above, the choice of a good range 
of values may not be trivial in particular for the original score. We investigate below in detail 
how the performance of TIGRESS depends on the scoring method and on these parameters R, 
a and L. 



°area 



with 

4>area(h) 



I. 



2.5 Other GRN inference methods 

We experimentally compare TIGRESS to several other GRN inference methods. We use the 
MATLAB implementations of CLR [ ] and GENIE3 [19]. We run ARACNE [26] using the R 
package minet. We keep default parameter values for each of these methods. Results borrowed 
from the DREAM5 challenge [23] were directly obtained by each participating team. 

2.6 Performance evaluation 

Given a gene expression data matrix, each GRN inference method outputs a ranked list of 
putative TF-TG regulations. Taking only the top K predictions in this list, we can compare 
them to known regulations to assess the number of true positives (TP, the number of known 
regulations in the top K predictions), false positives (FP, the number of predicted regulations 
in the top K which are not known regulations), false negatives (FN, the number of known 
interactions which are not in the top K predictions) and true negatives (TN, the number of 
pairs not in the top K predictions which are not known regulations). We then compute classical 
statistics to summarize these numbers for a given K, including precision (TP/ (TP + FP)), 
recall (TP / (TP + FN)), and fall-out (FP / (FP + TN)). We assess globally how these statistics 
vary with K by plotting the receiver operating characteristic (ROC) curve (recall as a function 
of fall-out) and the precision-recall curve (precision as a function of recall), and computing the 
area under these curves (respectively AUROC and AUPR) normalized between and 1. 

For the datasets of DREAM5, we further compute a P- value for the AUROC and AUPR 
scores, based on all DREAM5 participants' predictions. This method, which was used by the 
DREAM5 organizers to rank the teams, involves randomly drawing edges from the teams' predic- 
tion lists and computing the probabilities of obtaining an equal or larger AUPR (resp. AUROC) 
by chance. More precisely, random lists are constructed as follows: for each row of the predicted 
list, an edge at the same position is drawn at random from all predictions. For an ensemble 
of such random lists, the areas under the curves are computed, allowing to estimate a random 
distribution. P-values were obtained by extrapolating the resulting histogram. We refer to [23] 
for more details on this scoring scheme. Finally, we compute the so-called overall score for a 
GRN inference method by integrating the AUROC and AUPR P-values as follows: 

overall score = - In (pauprPAUROc) ■ (5) 

3 Data 

We evaluate the performance of TIGRESS and other GRN inference methods on four benchmark 
datasets, each consisting of a compendium of gene expression data, a list of known TF, and a 
gold standard set of verified TF-TG regulations which we ideally would like to recover from the 
expression data only. Expression data are either simulated or experimentally measured under a 
wide range of genetic, drug and environmental perturbations. Table 1 summarizes the statistics 
of these four networks. 

The first three benchmarks are taken from the DREAM5 challenge [23]. Network 1 is a 
simulated dataset. Its topology and dynamics were modelled according to known GRN, and the 
expression data were simulated using the GeneNetWeaver software [32]. We refer the interested 
reader to [25, 24] for more information on this network. The second and third benchmarks 
are Network 3 and Network 4 of the DREAM5 competition, corresponding respectively to real 
expression data for E. coli and S. cerevisiae. Note that we do not use in our experiments Network 
2 of DREAM5, because no verified TF-TG is provided for this dataset consisting in expression 
data for S. aureus. 



7 



Additionally, we run experiments on the E. coli dataset from [ ], which has been widely 
used as a benchmark in GRN inference literature. The expression data was downloaded from 
the Many Microbe Microarrays (M 3D ) database [ 2] (version 4 build 6). It consists in 907 
experiments and 4297 genes. We obtained the gold standard data from RegulonDB [J 5] (version 
7.2, May 6th, 2011) that contains 3812 verified interactions among 1525 of the genes present in 
the microarrays experiments. 

As a pre-processing step, we simply mean-center and scale to unit variance the expression 
levels of each gene within each compendium. 



Network 


tt TF 


jj Genes 


Jj Chips 


t) Verified interactions 


DREAM5 Network 1 (in-silico) 


195 


1643 


805 


4012 


DREAM5 Network 3 (E. coli) 


334 


4511 


805 


2066 


DREAM5 Network 4 (S. cerevisiae) 


333 


5950 


536 


3940 


E. coli Network from [11] 


180 


1525 


907 


3812 



Table 1: Datasets used in our experiments. 



4 Results 

4.1 DREAM5 challenge results 

In 2010 we participated to the DREAM5 Network Inference Challenge, an open competition to 
assess the performance of GRN methods [23]. Participants were asked to submit a ranked list 
of predicted interactions from four matrices of gene expression data. At the time of submission, 
no further information was available to participants (besides the list of TF), in particular the 
"true" network of verified interactions for each dataset was not given. After submissions were 
closed, the organizers of the challenge announced that one network (Network 1) was a simulated 
network with simulated expression data, while the other expression datasets were real expression 
data collected for E. coli (Network 3) and S. cerevisiae (Network 4), respectively. Teams were 
ranked for each network by decreasing overall score (5), and an overall ranking was proposed 
based on the average of the overall scores over the three networks. 

We submitted predictions for all networks with a version of TIGRESS which we did not 
try to optimize, which we refer to as Naive TIGRESS below. Naive TIGRESS is the variant of 
TIGRESS which scores candidate interactions with the original score (3) and uses the arbitrarily 
fixed parameters a = 0.2, L = 5, R\ = 4, 000, R3 = R4 = 1, 000, where Ri refers to the number 
of runs for network i. The number of runs were simply set to ensure that TIGRESS would finish 
within 2 days on a single-core laptop computer. R\ is larger than R3 and R4 because the size of 
network 1 is smaller than that of networks 3 and 4, implying that each TIGRESS run is faster. 
The choice a = 0.2 followed previous suggestions for the use of stability selection [29], while the 
choice L = 5 roughly corresponded to the largest value for which no TF-TG pair had a score of 
1. 

Naive TIGRESS was among the top GRN prediction methods at DREAM5, ranking second 
among 29 participating teams in the in silico network challenge, and third overall. Table 2 
summarizes the results of the first three teams in average overall score. The winning method, 
both m silico and overall, was the GENIE3 method of [19]. GENIE3 already won the DREAM4 
challenge, confirming its overall state-of-the-art performance. It had particularly strong perfor- 
mance on the in silico network, and more modest performance on both in vivo networks. The 
ANOVA-based method of ['~ ] ranked second overall, with particularly strong performance on 
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Teams 


Network 1 


Network 3 


Network 4 




AUPR 


AUROC 


Score 


AUPR 


AUROC 


Score 


AUPR 


AUROC 


Score 


GENIE3 [19] 


0.291 


0.815 


104.65 


0.093 


0.617 


14.79 


0.021 


0.518 


1.39 


ANOVA-based [20] 


0.245 


0.780 


53.98 


0.119 


0.671 


45.88 


0.022 


0.519 


2.21 


Naive TIGRESS 


0.301 


0.782 


87.80 


0.069 


0.595 


4.41 


0.020 


0.517 


1.08 



Table 2: AUPR, AUROC and minus the logarithm of related p- values for all DREAM5 Networks 
and the three best teams. 



the E. coli network. Naive TIGRESS ranked third overall, with particularly strong performance 
on the in silico network, improving over GENIE3 in terms of AUPR. 

Interestingly, GENIE3 and TIGRESS follow a similar formulation of GRN inference as a 
collection of feature selection problems for each target gene, and use similar randomization- 
based techniques to score the evidence of a candidate TF-TG regulation. The main difference 
between the two methods is that GENIE3 aggregates the features selected by decision trees, 
while TIGRESS aggregates the features selected by LARS. The overall good results obtained by 
both methods suggest that this formalism is particularly relevant for GRN inference. 

4.2 Influence of TIGRESS parameters 

In this section, we provide more details about the influence of the various parameters of TI- 
GRESS on its performance, taking DREAM5 in silico network as benchmark dataset. Obvi- 
ously the improvements we report below would require confirmation on new datasets not used to 
optimize the parameters, but they shed light on the further potential of TIGRESS and similar 
regression-based method when parameters are precisely tuned. 

Starting from the parameters used in Naive TIGRESS {R = 4, 000, a = 0.2 and L = 5, 
original score), we assess the influence of the different parameters by systematically testing the 
following combinations: 

• original (3) or area (4) scoring method; 

• randomization parameter a G {0, 0.1 . . . , 1}; 

• length of the LARS path L G {1, 2 ... 20}; 

• number of randomization runs R 6 {1, 000; 4, 000; 10, 000}. 

Figure 2 summarizes the overall score (5) obtained by each combination of parameters on Net- 
work 1. 

A first observation is that the area scoring method consistently outperforms the original 
scoring method, for any choice of a and L. This suggests that, by default, the newly proposed 
area score should be preferred to the classical original score. We also note that the performance 
of the area score is less sensitive to the value of a or L than that of the original score. For 
example, any value of a between 0.2 and 0.8, and any L less than 10 leads to an overall score of 
at least 90 for the area score, but it can go down to 60 for the original score. This is a second 
argument in favor of the area scoring setting: as it is not very sensitive to the choice of the 
parameters, one may practically more easily tune it for optimal performance. On the contrary, 
the window of (a, L) values leading to the best performance is more narrow with the original 
scoring method, and therefore more difficult to find a priori. The recommendation of [29] to 
choose a in the range [0.2, 0.8] is clearly not precise enough for GRN inference. The best overall 
performance is obtained with (a = 0.4, L = 2) in both scoring settings. 
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Area (1 ,000 runs) 



Original (1,000 runs) 




0.3 0.5 0.7 
Area (10,000 runs) 



0.1 0.3 0.5 0.7 0.9 
Original (10,000 runs) 




Figure 2: Overall score for Network 1. From top to bottom, plots show the results for R = 1, 000, 
R = 4, 000 and R = 10, 000 for both the area (left) and the original (right) scoring settings, as 
a function of a and L. 



Regarding the relationship between a and L, we observe in Figure 2 a slight positive cor- 
relation for the optimal L as a function of a, particularly for the area score. For example, for 
R = 10 4 , L = 2 is optimal for a < 0.4, but L > 4 is optimal for a > 0.8. The effect is even more 
pronounced for R = 4, 000. This can be explained by the fact that when a increases, we decrease 
the variations in the the different runs of LARS and therefore reduce the diversity of features 
selected; increasing the number of LARS is a way to compensate this effect by increasing the 
number of features selected at each run. Another way to observe the need to ensure a sufficient 
diversity is to observe how the best parameters L and a vary as a function of R (Figure 3). It 
appears clearly that the optimal number of steps L* decreases when the number of resampling 
runs increases and stabilizes at 2. This is not a surprising result. Indeed, when more resampling 



10 



is performed, the chance of selecting a given feature increases. The number N of non zero scores 
subsequently increases and it thus becomes unnecessary to look further in the regularization 
path. On the other hand, the value of a* lies steadily between 0.3 and 0.5, suggesting that the 
adjustment to the number of bootstraps can mostly be made through the choice of L. 
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Figure 3: Optimal values of parameters L, a and N with respect to the number of resampling 
runs 

Finally, we unsurprisingly observe that increasing the number R of resampling runs leads 
to better performances. On Figure 4, we show the overall score as a function of R with L = 2 
and a = 0.4. We clearly see that, for both scoring methods, increasing the number of runs is 
beneficial. The performance seems to reach an asymptote only when R becomes larger than 
5,000. 

4.3 Comparison with other methods 

Figure 5 depicts both the ROC and the Precision/Recall curves for several methods on Network 
1. Table 3 summarizes these performances in terms of AUPR, AUROC and related p-values as 
well as the overall score. Here, TIGRESS was run with a = 0.4, L = 2 and R = 8, 000 which 
corresponds to the best performance of the algorithm, as investigated in the previous section. 

TIGRESS outperforms all methods in terms of AUPR and all methods but GENIE in terms 
of AUROC. Moreover, the shape of the Precision/Recall curve suggests that the top of the 
prediction list provided by TIGRESS contains more true edges than other methods. The ROC 
curve, on the other hand, focuses on the entire list of results. Therefore, we would argue that 
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Figure 4: Overall score as a function of R. In both scoring settings, a and L were set to 0.4 and 
2, respectively. 
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Figure 5: ROC (left) and Precision/Recall (right) curves for several methods on Network 1. 
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Method 


AUPR 


PAUPR 


AUROC 


PAUROC 


Overall score 


TIGRESS 


0.320 


1.17e-145 


0.789 


3.74e-67 


105.68 


GENIE3 


0.291 


1.60e-104 


0.815 


3.06e-106 


104.65 


Naive TIGRESS 


0.301 


7.20e-118 


0.782 


3.48e-59 


87.80 


CLR 


0.265 


1.82e-73 


0.782 


1.41e-58 


65.30 


ANOVA-based 


0.245 


8.17e-53 


0.780 


1.34e-56 


53.98 


ARACNE 


0.276 


1.73e-85 


0.672 


9.82e-01 


42.38 



Table 3: Comparison of different methods on Network 1 of the DREAM5 challenge. The perfor- 
mance of GENIE3, Naive TIGRESS and ANOVA were obtained during the DREAM5 competi- 
tion. TIGRESS corresponds to the choice of parameters leading to the best performance (area 
score, a = 0.4, L = 2, R = 8,000). We ran CLR and ARACNE using public implementations 
of these methods. 

TIGRESS is more reliable than GENIE in its first predictions but contains overall more errors 
when we go further in the list. 

4.4 In vivo networks results 

Since Naive TIGRESS did not perform very well on the in vivo networks at the DREAM5 
competition (Table 2), we now test on these networks TIGRESS with the best parameters 
selected on the in silico (area score, a = 0.4, L = 2 and R = 10,000). Table 4 shows the values 
of AUPR, AUROC, related p- values and overall score for DREAM5 networks 3 and 4 reached 
by TIGRESS. 



Network 


AUPR 


PAUPR 


AUROC 


PAUROC 


Overall score 


DREAM5 Network 3 


0.0660 


4.79e-06 


0.5887 


6.66e-02 


3.25 


DREAM5 Network 4 


0.0199 


5.86e-01 


0.5143 


2.02e-01 


0.46 



Table 4: TIGRESS performance on DREAM5 Networks 3 and 4. 



The results on these two networks are overall disappointing: TIGRESS does not do better 
than Naive TIGRESS. In fact, both sets of results are very weak. Without attempting to re- 
optimize all parameters for each network, one may wonder whether the parameters chosen using 
the in silico network are optimal for the in vitro networks. As a partial answer, Figure 6 shows 
the behavior of the overall score with respect to L for Networks 3 and 4. Interestingly, it seems 
that a larger L is preferable in this case, suggesting that one may have to adapt the parameters 
to the size of the network. Indeed, networks 3 and 4 contain respectively 4, 511 and 5, 950 nodes, 
making them much larger than the in silico network we tuned the parameters on. However, the 
improvement is not dramatic in absolute value. 

On Figure 7 we compare Precision/Recall and ROC curves obtained with TIGRESS with 
several other algorithms on the E. coli network from [11]. Table 5 compares the areas under the 
curves. TIGRESS is comparable to CLR, while GENIE3 outperforms other methods. However 
the overall performance of all methods remains disappointing. 

4.5 Analysis of errors on E. coli 

To understand further the advantages and limitations of TIGRESS, we analyse the type of 
errors it typically makes taking the E. coli dataset as example. We analyse FP, i.e. cases where 
TIGRESS predicts an interaction that does not appear in the gold standard GRN. 
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Figure 6: Overall score with respect to L for networks 3 and 4 (a = 0.4, R = 10, 000). 




Figure 7: Precision/Recall (Left) and ROC (Right) curves for several methods on the E. coli 
dataset. 



Method 


AUPR 


AUROC 


TIGRESS 


0.0624 


0.6026 


ARACNE 


0.0498 


0.5531 


CLR 


0.0641 


0.6019 


GENIE3 


0.0814 


0.6375 



Table 5: TIGRESS compared to state-of-the-art methods on the E. coli Network. 



We focus in particular on quantifying how far a wrongly predicted interaction is from a true 
one, and introduce for that purpose the notion of distance between two genes as the shortest 
path distance between them on the gold standard GRN, forgetting about the direction of edges. 
For two genes Gl and G2, we call G1-G2 a distance-x link if the shortest path between Gl 
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and G2 on the true network has length x. Figure 8, shows the distribution of these distances 
for spuriously discovered edges over the gold standard network, i.e. the actual proportion of 
distance- a; links, with x 6 {1,2,3,4, > 4}. We write p x the proportion of spurious TG-TF 
couples with distance x. 




1 2 3 4 >4 
Length of shortest path 



Figure 8: Exact distribution of the shortest path between spuriously predicted TF-TG couples. 

Figure 9 depicts the distribution of distance- x proportions among the spuriously detected 
edges, as a function of the number of predicted edges. Dotted lines represent the 95% confidence 
interval around the exact distribution (p x ) x . For a given number r of spuriously predicted edges, 
this interval is computed as 

Q0.025(Px) _ 90.975 fe) 
r r 

where q a {Px) represents the quantile of order a of a hypergeometric distribution %{Ns,p x Nsir) 
and N$ is the total number of spuriously predicted edges. 

We observe that most of the recovered false positives appear as distance-2 edges in a signif- 
icantly higher proportion than p2 whereas significantly less distance-> 4 edges are discovered. 
These results strongly suggest that most of TIGRESS errors - especially at the top of the list 
- are indeed sensible guesses, where the two nodes, spuriously discovered with a parent/child 
relationship are actually separated by only one other node. In Table 6, we detail the three 
possible patterns observable in this situation. 

Figure 10 focuses on distance-2 errors. Note that some edges show more than one pattern, 
e.g. the first spurious edges are both siblings and couples. It appears that most of them are 
siblings and can thus be interpreted as spurious feed-forward loops. We believe that this can 
be explained by three main reasons: i) the discovered edges actually exist but have not been 
experimentally validated yet; ii) there is more of a linear relationship between siblings than 
between parent and child; iii) some nodes have very correlated expression levels, making it 
difficult for TIGRESS to tell between the parent and the child. 
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Figure 9: Distribution of the shortest path length between nodes of spuriously detected edges and 
95% confidence interval for the null distribution. These edges are ranked by order of discovery. 



16 



Name 


Illustration 


Description 


Siblings 


A 


Gl and G2 have a common 
parent. They are siblings. 


Couple 


© © 

o 


Gl and G2 have a common 
child. They are a couple. 


Grandparent / Grandchild 


(H>-O0 


Gl has a child that is a parent 
of G2. Gl is a grandparent of 
G2. 



Table 6: Distance-2 patterns between two nodes Gl and G2 in a directed graph. 




Figure 10: Distribution of distance 2 errors. 95% error bars were computed using the quantiles 
of a hypergeometric distribution. 
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5 Discussion 



In this paper, we presented TIGRESS, a new method for GRN inference. TIGRESS expresses 
the GRN inference problem as a feature selection problem, and solves it with the popular LARS 
feature selection method combined with stability selection. It ranked in the top 3 GRN inference 
methods at the 2010 DREAM5 challenge, without any parameter tuning. We clarified in this 
paper the influence of each parameter, and showed that further improvement may result from 
finer parameter tuning. 

We proposed in particular a new scoring method for stability selection, based on the area 
under the stability curve. It differs from the original formulation of [29] which does not take into 
account the full distribution of ranks of a feature in the randomized feature selection procedure. 
Comparing the two, we observed that the new area scoring technique yields better results and 
is less sensitive to the values of the parameters: practically all values of, e.g., the randomization 
parameter a yield the same performance. Similarly, the choice of the number L of LARS steps 
to run seems to have much less impact on the performance in this new setting. As we showed, 
the original and area scores for a feature t can be both expressed in a common formalism as 
E [<j)(H)\ for different functions <j), where Ht is the rank of feature t as selected by the randomized 
LARS. It could be interesting to systematically investigate variants of these scores with more 
general non-increasing functions (j), not only for GRN inference but also more generally as a 
generic feature selection procedure. 

Comparing TIGRESS - as tuned optimally - to state-of-the-art algorithms on the in silico 
network, we observed that it achieves a similar performance to that of GENIE3 [19], the best 
performer at the DREAM5 challenge. However, TIGRESS does not do as good as this algorithm 
on in vivo networks. GENIE3 is also an ensemble algorithm but differs from TIGRESS in that 
it uses a non-linear tree-based method for feature selection, while TIGRESS uses LARS. The 
difference in performance could be explained by the fact that the linear relationship between 
TGs and TFs assumed by TIGRESS is far-fetched given the obvious complexity of the problem. 

A further analysis of our results on the E. coli network from [J 1 ] showed that many spuriously 
detected edges follow the same pattern: TIGRESS discovers edges when in reality the two nodes 
are siblings, and thus tends to wrongly predict feed-forward loops. This result suggests many 
directions for future work. Among them, we believe that operons, i.e. groups of TGs regulated 
together could be part of the problem. Moreover, it could be that there is more of a linear 
relationship between siblings than between parent and child, as TFs are known to be operating 
as switches, i.e. it is only after a certain amount change in expression of the TF that related 
TGs are affected. However, it is worth noting that in vivo networks gold standards may not be 
complete. Therefore, the hypothesis that TIGRESS is actually correct when predicting these 
loops cannot be discarded. 

While it seems indeed more realistic not to restrict underlying models to linear ones, it is fair 
to say that no method performs very well in absolute values on in vivo networks. For example, 
performances on the E. coli network seem to level out at some 64% AUROC and 8% AUPR 
which cannot be considered satisfying. This suggests that while regression-based procedures such 
as TIGRESS or GENIE3 are state-of-the-art for GRN inference, their performances seem to hit 
a limit which probably cannot be outdistanced without some changes in the global approach 
such as adding some supervision in the learning process as, e.g., investigated in [30]. 
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